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Cosmic ray propagators 
in the fractional differential model of bounded anomalous diffusion 
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Summary. Fractional differential approach to cosmic ray physics problems is discussed. A short review in this field is given, some 
results are represented, analyzed and criticized. A new model called the bounded anomalous diffusion model is offered. Its equation 
includes the fractional material derivative which allows to take into account the finite speed of cosmic rays particles. 

Introduction 

For a long time, cosmic rays (CR) propagation through the Galaxy have been described in frame of the normal diffusion 
theory, based on the equation 

§- t ~CA ^(r,t) = 5(t)5(r), 

where C is the diffusion coefficient independent on frequency. This equation is derived under the assumption that fluc- 
tuations of interstellar magnetic fields, giving particle trajectories a specific (Brownian) character, are characterized by 
definite space sizes and look more homogeneous at large scales. 

Rising attention to fractional calculus and its applications haven't passed CR physics as well (see [Uchaikin(2008)|). 
An operator only symbolically different from the fractional derivative with respect to time, appeared (being unrecog- 



nized) in fChuvilgin & Ptuskin (1993)j. The same derivative we meet in the equations derived by |Ragot & Kirk (1997) 



Balescu (1995)] (see also the reviewing part of recent work [Webb et al. (2006) |). Replacing the first time-derivative 
in the diffusion equation by its fractional analog reflects slowing-down of CR diffusion caused by magnetic traps (see 
[Potman (1975)| ). This kind of diffusion process was called subdiffusion. 

At the same time, the turbulent (fractal) character of the interstellar magnetic field (see [Ruzmai kin et al. (1988)| ?]) 
accelerates the diffusion. The turbulent diffusion {superdiffusion) was described by [Moni n~("l955)| by means of diffu- 
sion equation with fractional (2/3) power of Laplacian. [Saichev & Zaslavsky (1997) | and | Uchaikin & Zolotarev (19 99) | 
combined both of these regimes (super- and subdiffusions) to one anomalous diffusion process, described by the space- 
time fractional equation 

Df + C(-A) Q / 2 J <A(r, t) = ^L-S{t), (1) 

Its solution depends on two parameters: fractional orders of the partial derivatives with respect to the coordinates (a 6 
(0, 2]) and time (j3 6 (0, 1]). The characteristic features of this solution are the presence of heavy power-law tails 
(for a < 2) and the law of diffusion packet spreading t^l a . The family of distributions was named fractional stable 
distributions by [Kolokoltsov et al. (2000)], and investigated later in detail. 

First application of fractional diffusion model to CR physics was connected with energy spectrum problem (Lagutin 
and Uchaikin, [Lagutin et al. (2001) Lagutin & Uchaikin (2001) Lagut in & Uchaikin (2003)) ). The fractional diffusion 
model was used because interstellar magnetic field heterogeneities were testified to take large-scaled (fractal) charac- 
ter (see [Kulakov & Rumjantsev (1994) |). The supernova remains analysis shows the presence in this region of gas com 



ponents with different physical parameters (T e ~ 5 -j- 10 K, n e ~ 0.1 -j- 10 m _d ), what can be the sequence of extreme 
heterogeneity of interstellar medium. These facts and other data concerning the heterogeneities of matter density p and 
magnetic field intensity H oc p q , q ~ 1/3 -i- 1/2 ([Ruzmaikin et al. (1988) |) within the length range 100-150 pc, giving 
rise to uncertainties of diffusion model to be applicable to CR transfer description, stimulated the use of the superdiffusion 
model, based on the fractional Laplace operator: 



! + c(£x-Ar/ 2 



${r,t,E) = S(r,t,E). 



(2) 



This model seemed to be able to explain the experimentally observed "knee" in the energy spectrum, i. e. increasing of 
the exponent 77 in power representation of the E~ v spectrum while passing from the region ~ 10 2 Gev/nucleon into that 
of ~ 10 5 Gev/nucleon. The explanation was connected with the presence of power kind asymptotics of the solution at 
small and large distances from a point source. Moreover, this was in accordance with the self-similarity hypotheses which 
led Monin to the same equation in the turbulent diffusion problem. 
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Later, Lagutin et al. published a number of works in which the distributions i/j(r,t) obtained by 
| Saichev & Zaslavsky (1997) Uchaiki n & Zolotarev (1999)) were compared with experimental data by choosing the pa- 



rameters. In the process of those investigations, the parameters were changed from initial values a = 1.7, /3 = 1 -j- 0.8 
used in our first works in 2001 to a = 0.3 and (3 = 0.8 in 2004 (see | Lag utin & Tyumentsev (2004)| , p. 13); thus, the 



ratio /3/a increased from 0.47 to 2.67. The former values were not unreasonable, whereas the latter values seem to be 
absurd. Indeed, in this case, the cloud of CRs instantaneously emitted by a point source spreads according to the law 



whereas the limit speed of the particles is the light speed c and the diffusion packet cannot spread faster then R c = ct. 
The cause of such conclusion is that Eqs. (1) and (2) describe processes of instantaneous jumps separated by random 
periods of immobility state. This contradiction to real physical process was first overlooked because of the absent of a 
clear interpretation of fractional derivatives. 

On fractional derivatives and fractals 

The pecularity of fractional derivatives is their separation from the usual differential and the related notion of increment. 
The representation that the increment Ay is for some reasons proportional to ((Ai)") does not concern the Riemann- 
Liouville derivative and its three-dimensional Riesz generalization. In this case, one has to consider a special form of the 
increment, the difference of the fractional (f-th) order: 



k=0 ^ ' 

However, this difference cannot be interpreted in an explicit way, when v is not an integer. For this reason, the derivation of 
fractional equations usually begins with integral relations, which are then subjected to the Fourier-Laplace transformations 
r i-> k, tn>A and asymptotic expansions. As a result, the products |k| Q /(k, A), a e (0, 2], and A^/(k, A), /3 S (0, 1], 
are treated as the Fourier transform of the fractional Laplacian and the Laplace transform of the Riemann-Liouville 
fractional derivative of the orders a and /?, respectively, 

oo 

|k|«/(k, ■) = J e ikr A«/ 2 /(r, -)dv, A^/0, A) = J Df /(•, t)dt. 

o 

The very expressions |k| Q /(k, •) and A' 3 /(-,A) appear as asymptotic expressions (for k —> 0,A — > 0)), and inverse 
transformations using Tauberian theorems lead to the power-law (in r and t) functions in the asymptotic limit of large 
values of arguments. Such distributions characterize, in particular, fractal structures and processes whose main property is 
self-similarity. For example, if a sequence of random points is located on the real axis so that the distances Rj between the 
neighboring points are identically distributed and independent and the distribution of points is self similar on average (i.e., 
(N(x)) = (N{l))x v , v 6 (0, 1]), the equation for the distribution density of Rj includes the fractional derivative of the 
order v. The resulting distribution of random points has all of the attributes of a stochastic fractal, including intermittence. 
As v increases, these attributes weaken and the equation with v = 1 is the first-order equation whose exponential solution 
provides the model of an inhomogeneous Poisson ensemble, which has a homogeneous form at large scales. This is the 
relation of fractal structures with a fractional derivative. 

Since the most important property of space plasma is turbulence, manifesting self-similar structures of a power-law 
(fractal) type, the application of the fractional derivative technique to diffusion in this medium seems to be appropriate in 
this case. 

Application of the fractional differential equation (1) to anomalous diffusion is based on the continuous time random 
walk model introduced by [Montroll & Weiss (1965)]. In this model, the walk of a particle is represented by a sequence 
of instantaneous jumps with random lengths at random times between which the particle is at rest. The lengths of the 
jumps and intervals during which the particle is at rest (in "a trap") are independent of each other. To illustrate the relation 
of the walk scheme with the real transport of CRs in the galactic magnetic field, let us divide space into cubic cells and 
specify each particle inside the i-th cell by the radius vector of the center of this cell (Fig. 1). At a random time T after 
entering this cell, the particle passes to one of the six neighboring cells and the corresponding vector jumps to the center 
of this new cell at the time of the intersection of the interface between the old and new cells. After a random time T' , the 
particle passes to the next neighboring cell and the vector again instantaneously jumps, etc. If these cells were identical in 
the properties, the walk of the specifying vector would be discretized Brownian motion, which is transferred to ordinary 
diffusion with an increase in the scale. In case of a regular medium all cells are identical and free motion regime is absent. 
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Figure 1 : The coarse-grained trajectory in regular (left panel) and fractal (right panel) medium. 



However, the strongly turbulent high magnetic field is not inherent in each cell. According to the current representation, 
which was formed five decades ago, most of the space between magnetic clouds is filled with lower quieter fields whose 
smooth field lines can hold individuality at a long distance. The charged particles of CRs move along these field lines 
in helical paths and from time to time enter trap clouds, where they can stay for a long time and "forget" their initial 
direction of motion. In this case, a jump from one trap to another is not instantaneous as in the above case: the particle 
rather intersects "almost empty" cells, covering large distances R. The distribution of R in a fractal medium involves 
a long power-law tail. These transitions require the consideration of a finite velocity of the particle, more precisely, the 
leading center of the particle. 

The bounded anomalous diffusion model 



The most important consequence of the finiteness of the velocity of free particles is the finiteness of the spatial distribution: 
the probability density beyond a sphere with the radius vt and the center at the instantaneous point source is zero in this 
case. Let us refer to this process as bounded anomalous diffusion in order to distinguish it from unbounded anomalous 
diffusion in which the motion of a particle is represented as a sequence of instantaneous jumps from one point of space to 
another: the particle arrives at the latter point at the same time at which it leaves the former point at any distance between 
these two points. The delay time (in a trap) is not related to this distance and to the motion as a whole. If the traps are 
removed from this model, it becomes senseless, because the particle instantaneously flies to infinity and leaves the system 
under consideration. The bounded anomalous diffusion model involves not only the time spent in the traps, but also the 
time taken for the motion of the particle and, for this reason, is meaningful even in the absence of traps. Finally, since the 
bounded anomalous diffusion propagator vanishes beyond the sphere with the radius vt, all of its moments are finite. 
The effect of the finiteness of the velocity of free motion on the walk process described above, which was 
investigated in ||Uchaikin (1998a)| |Uchaikiii t (1998b)] |Uchaikin (1998c)] |Uchaikin (1998d)] |Zolotorev et al. (1999)| 



|Uchaikin & Yarovikova(2003)| , significantly changes the continuous time random walk model. In this case, the parti- 
cle at the observation time can be in one of two states, rest and motion. The corresponding components of the probability 
density are denoted as ipo (r, t) and ipi ( v , t) so that the total density is 

^(r,t)=^ (r 1 i)+^ 1 (r,t). 

The rates of the 1 — >• and — > 1 transitions per unit volume near the point r are denoted as Fi_ > o(r, t) and i*o^i( r ; t), 
respectively. It is obvious that the particle passing to the state of rest at the point r at the time t — t' remains in this state at 
the observation time t with the probability Q(t') = L q(t)dt, and the particle leaving the trap at the point rr' intersects 
a unit area at the point r without an interaction with the probability -P(r') = J* °° p(r' + £fi)d£, fl = r'/r'. Since this 
transition takes r'/v seconds, 

oo 

#•,*) = J dtQWF^ofrt-^ + v- 1 J dr'Pir^Fo^ir-r^t-r'/v), (3) 
o 

The factor v^ 1 appears in front of the integral because the integral gives the flux of particles, whereas ipi is the concen- 
tration of these particles. The transition rates (if the particle begins its evolution in a trap at the origin of the coordinates 
at the initial time) are related as 

Fi_x,(r,t)= j dr'p{r')F ^ 1 {Y-T l ,t-r'/v) + 5{r)5{t), (4) 



4 



F _n(r,i) = J dTq{T)F^ (T,t-T), 



The FourierLaplace transformation 

oo 

ip(r,t) ^ V(k, A) = e ikr - A V(r,*)rfi 

o 

reduces the system of Eqs. (3)-(5) to the form 

L„(A,k)V>(k,A) = [1 -p(k,A/ W )g(A)]V(k,A) = S„(A,k), 

where 



(5) 



and 



S„(A,k) = 0(A) + (l/v)P(k, X/v)q(X) 



p(k,A/u) - y p(r)e- (A/u)r e ikr dr. 



The quantity P(k, A/w) is determined similarly. If the tails of the distributions of R and T are of a power-law character 
with exponents a 6 (0, 2) and /3 e (0, 1], respectively, 



P(R > r) 



-r , r — > oo, 



P(T > t) = Q(t) 



r(l-a) 

Then, using Tauberian theorems, one can show that for k — > 0, A — > 



r(i-/3) 



t -> oo. 



L„(A,k) ~BA^ + 



(i?)(A/w-ikfi)+A((A/w-ikn) a ), a>l, 
A((A/w -ik$7) a ), a<l. 



where f2 is a random direction of the motion, which is assumed to be isotropically distributed. The comparison of different 

terms at A — > provides the following conclusions. 

The asymptotic expressions at v = oo for both cases have the same form 



L„(A,k) ~B[A^ + C 00 |k|' 



Coo =i4|cos(a7r/2)|/[(a + l)B], 



This expression provides fractional differential equation (1) of unbounded anomalous diffusion. At v < oo, a = 2, (3 = 1 
(in view of the mentioned isotropy, (ft) = 0), the transport operator has the form 

L v (\, k) ~ BX 13 + (R)(X/v - ikfl) + A((X/v - ikfl) 2 ) = (A/v 2 )X 2 + (B + (R)/v)X - (A/3)k 2 , 

corresponding to the telegraph equation 



A_cP_ 

v 2 dt 2 



B 



A A 
v J dt 3 



^(r,t) = S v (r,t), 



which describes the bounded normal diffusion. Under the same conditions at < 1, the following fractional variant of 
the subdiffusion telegraph equation is obtained: 



AcP_ 
^di 2 



(R)_d_ 
v dt 



i>(r,t) = S v (r,t). 



In both bases, the form of the "deep" time asymptotic expression for ^(r, t) is independent of the velocity and leads to the 
equations of bounded diffusion and subdiffusion, respectively. The same property is observed at a < 1 when f3 < a. If 
a < f3 < 1 (this relation between the parameters is used by Lagutin and Tyumentsev), the asymptotic equation for ip(r, t) 
has the form unusual for a diffusion process with a pseudodifferential operator (D t + v fiV) Q averaged over the directions 



(A/v a )((D t + vflV) a )i>(v,t) =S v (r,t). 



Here, the operator 



(D t + w nv) v v(r, t) = T{1 1 _ [/) (Dt + ««V) J 



tp(r — vft(t — t),t) 



dr 
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is the fractional generalization of the material derivative that is agree with results obtained by [Sokolov & Metzler(2003)| 
for one-dimensional Levy walks. 

If we consider the case of homogeneous distribution of particles ip(r, t) = ip(t) this operator becomes the fractional 
Rieman-Liouville derivative with respect to time: 

(D t + vSlV)"<tp(r,t) oD^(t) 

In the stationary problem, when the function tp(r, t) = -0(r) does not depend on time, 

(Dt + vilV)" i/j(r) i-> («nV)^(r), 

this operator represents fractional generalization of directional derivative. 

The dependence on the velocity remains at any time, but the dependence on disappears. The distribution ip(r,t) cor- 
responding to this case has a specific U shape in the region bounded by the radius vt beyond which it vanishes (see 
QZolotorev et al. (1999)||Uchaikin & Yarovikova(2003)) ). 

The bounded anomalous diffusion propagator 

In one-dimensional case, for a < 1, the equation of bounded anomalous diffusion without traps takes the form: 



dt 



d_ 

dx 
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dt V dx 



r 



2T(1 - v) 



[S(x - vt) + 6{x + vt)}, 



Solutions of this equation can be expressed through elementary functions (see [Uchaikin & Sibatov(2009)|) 



2 sin ttv 



(l-x 2 /v 2 t 2 ) 1 



(1 - x/vt) 2v + (1 + xjvt) 2v + 2(1- x 2 /v 2 t 2 ) v COS TXV ' 



Solutions can also be written in terms of fractional stable densities (see | Uchaikin & Sibatov(2009)|), that are useful for 
probabilistic interpretation of these distributions. 

In Fig. 2 (left panel), the analytical solutions (lines) are compared with the results of Monte Carlo simulated random 
walks with finite velocity of motion (points). In Fig. 2 (right panel), the influence of ballistic restriction on the dynamics 
of diffusion packet spreading is demonstrated schematically. 

To make clear the role of correlations between path lengths and waiting times in the model under consideration, we 
compare propagators of bounded and unbounded anomalous diffusion. In the latter model we take identical distributions 
for waiting times and path lengths (a — (3). In the model of unbounded anomalous diffusion path lengths and waiting 
times in traps are independent. We take the model of bounded anomalous diffusion without traps, time delay is provided 
by finiteness of rays propagation velocity (v — 1). The results for the one-dimensional case are presented in Fig. 3 
(left panel), and for the 3D-model in the right panel. From the graphs, we can see that distinction in kind between 
propagators disappears only when a = f3 > 2. For values a — (3 < 1, distinctions are very strong, shapes of packets, 
spreading laws, behaviors near the ballistic boundaries and near the source are essentially different. When 1 < a = 
(3 < 2, distinctions are also sufficient despite the fact that mean path length is finite. In one case distributions are 
bounded, in another case they are unbounded. In the model of bounded anomalous diffusion, the front near |r| = vt 
appears. The densities differ quantitatively near the source as well. All these facts say that the results obtained on the 



base of unbounded anomalous diffusion and presented in the works |Lagutin et al. (2001) Lagutin & Uchaikin (2001) 



Lagutin & Uchaikin (2003) Lagutin & Tyumentsev (2004) Lagutin et al. (2009) | need to be re-examined 



Concluding remarks 



The work (| |Zolotorev et al. (1999)| , p. 1424) on anomalous diffusion ended with the phrase: "The last property can be 
a reason for the conclusion that the superdiffusion equation is inapplicable to the description of real physical processes 
with the characteristic exponents a < 1" . The cautious word "can" is not accidental here: if the walk of the particle 
is considered in the space of, e.g., velocities or momenta, an instantaneous jump of the particle at large "distance" is 
allowable and even justified in the framework of the commonly accepted model of instantaneous collisions. The use of 
the model with an infinite velocity of particles without any boundary conditions to describe the walk of the particles in the 
coordinate space at a < j3 makes both the results and conclusions doubtful (in particular, the conclusion that "the self- 
consistent description of the experimental data on the spectra of individual groups of nuclei and the spectrum of all of the 
particles is not achieve" [Lagutin et al. (2009) | p. 601). We hope that taking into account the finiteness of the propagation 



velocity of CRs in the framework of the bounded anomalous diffusion proposed in this work will be more fruitful. 
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Figure 2: Comparison of analytical solutions (lines) with results of numerical Monte Carlo simulation (points) in one-dimensional case 
(left panel). Evolution of the Gaussian diffusion packet and the propagator of the bounded anomalous diffusion model, a — 0.5 (right 
panel). Dashed lines represent the time dependencies of the diffusion packet width in both cases. 




\\i(x,t) ] 



0.01 



1E-3 



1E-4 



ot=p=1.50 



-100 



% 

% 

«> 



\|/ 3 (r,0 



1E-2 - 



1E-3 - 



♦oV°„ 

♦ ° 



a=p=1.5 



100 



100 



200 



Figure 3: Comparison of propagators of the bounded anomalous diffusion (circles) and unbounded diffusion (rhombs) for a = 0.5 and 
1.5 in the ID case (left panel) and 3D case (right panel). 
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